On multisolitonic decay behavior of internal gravity waves 
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Abstract 

We claim that changes of scales and fine-structure could increase from multisoliton behavior of 
internal waves dynamics and, further, in the so-called "wave mixing". We consider initial-boundary 
problems for Euler equations with a stratified background state that is valid for internal water waves. 
The solution of the problem we search in the waveguide mode representation for a current function. 
The orthogonal eigenfunctions describe a vertical shape of the internal wave modes and satisfy a 
Sturm-Liouville problem for the vertical variable. The Cauchy problem is solved for initial conditions 
with realistic geometry and magnitude. We choose the geometry and dimensions of the McEwan 
experiment with the stratification of constant buoyancy frequency. The horizontal profile is defined 
by numerical solutions of a coupled Korteweg-de Vries system. The numerical scheme is proved to be 
convergent, stable and tested by means of explicit solutions for integrable case of the system. Together 
with the solution of the Sturm-Liouville problem it gives the complete internal waves profile. 
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Introduction 



There are papers (see S. B. Leble (1991) for a review) in which the problem of nonlinear waves evolution 
in stratified medium is in some way connected with soliton theory, see, e.g., Novikov et al (1984). We 
mean the conditions in which the physical (or environmental) conditions need account of nonlinearity and 
dispersion which joint action is of essential (sometimes - inevitable) importance. Evolution of internal 
gravity waves (IGW) in geophysical hydrodynamics and its further decay to turbulent spots is important 
example of such phenomena, Pedlosky (1979). The scenarios of this transition intensely discussed in 
literature of atmosphere and ocean physics Gavrilov, Fukao, Nakamura, Tsuda (1996). There are also 
publications about laboratory experiments devoted to IGW dynamics in tanks Segur, Finkel (1985). The 
problem have a link to a general problems of hydrodynamics as illustrated by Zeitunyan (1999). It is the 
Euler equations used as a basis of a description in which modes with different scales are generated within 
the process of a wave propagation over an initial density stratification taking into account its possible 
dynamic changes. This question relates to the behavior of a solution within the limiting procedure when a 
dispersion coefficient go to zero. In our case of internal waves it means the zero buoyancy frequency limit 
(N — > ) that in turn is proportional to the vertical density gradient. Let us mention also investigations of 
similar problems with a small dispersion parameter initiated by Lax and Levermore (1983) for the model 
of a single Korteweg-de Vries (KdV) equation. The dispersion in such model is described by the term with 
third derivative by a coordinate (x). 



In this paper we study an evolution of a disturbance from an initial condition that model an excitation 
of the internal wave in the geometry and scales of McEwan (1983 a,b) experiments. We believe that the fine 
structures that appear in density distribution may be explained by a model based on a system of coupled 
Koreweg-de Vries (cKdV) equations. This system arises for coefficient functions related to a functional 
basis of transversal coordinate. It demonstrates the behaviors of solutions that have many features of 
multi-solitonic decay (see, again, Novikov et al (1984)) The wavetrains contain chains of such fluctuations 
of basic wave variables that propagate with velocity up to the most quick one. Our investigations show: 
interaction between the basic transverse guide modes does not destroy the general picture; only some phase 
changes appear. Of course the system of cKdV has higher order in time derivatives, the opposite directed 
wave train is launched, but the wave also has the multisoliton form. 

The derivation of the system of cKdV may be considered as the direct application of Galerkin - Bubnov 
projection procedures that is successfully used in numerical simulation of initial-boundary problems with 
a good mathematical justification (the basis of this method are in Gottlieb and Orszag (1989), see, also 
Moser (1966)). After the variables division in the linearized problem some spectral problem appears. Its 
eigen functions of transverse variables form a basis. Hence the (infinite) 1+1 cKdV system is equivalent 
to the 2+1 basic Euler system under consideration within some class of initial conditions. While numeric 
simulation delivering we cut this infinite system restricting ourselves by some mode number, comparing 
few such examples. We, however believe that further extending of the system does not change the general 
behavior of the wavetrains. We plan checking of the hypothesis in the forthcoming paper but more 
concentrate now on the cKdV model that is also interesting by itself (e.g. Hirota R., Satsuma J. (1981) , 
A. Perelomova (2000)). It was studied mathematically via many approaches as described recently by Ayse 
Karasu (1997). There are integrable examples (R. Dodd, A.Fordy (1982)): it helps to test the numerical 
scheme to be introduce in this article. 

For such systems we construct a finite-difference scheme (Sec. 3) and use it for numerical integration. 
The theorems about convergence and stability we formulate and prove here confirm the statements we 
made. The computations description we place in the section 4 and the proofs in Appendix. 

In this article we do not touch the mixing problems in its full extent as well as interaction with 
the stratification mode (background or stationary state) itself. Such problems solution need account of 
boundary conditions and simulations for a longer time intervals. Account of such ingredients in the model 
could shade the main phenomena we intend to demonstrate. Therefore we concentrate attention of a 
reader on the initial stage of the wave field development. Already within this initial stage period the 
horizontal multi-solitonic fine structure prevails: it seems very important. Our general approach and the 
numerical scheme power allow to consider longer time to account reflections phenomena. A comparison 
of the evaluated stream function fields with ones of the McEwan experiments show many similar features, 
that, seems to verify the model from the point of the real fluid dynamics. 

2 Formulation of the problem 

2.1 Laboratory experiment of McEwan 

Experiment were made using a rectangular plate-glass tank with dimensions shown in figure (1) below. 
This tank is filled by a linearly stratified salt solution. A wave-forcing dismountable paddle swung about 
horizontal lateral axis ia used ||. 



Axis of rotation 




50 cm 

Tank filled with a continuous stratified ambient state water 
Fig. (1) A. D. McEwan experimental configuration. 

Bring the fluid to a state of incipient breaking and photographing the evolution of these breaking events 
induced by internal wave action. McEwan gives his observation along the following stages: 

(a) Overturning. 

(b) Development of interleaving microstructure. 

(c) Static stability is restored but microstructure is preserved. 

(d) Gravitation to an equilibrium has changed the surrounding density profile between extremum isopyc- 
nals. 

We simulate only tile stage (b) above to show the obvious plots for the mixing phenomena. 



3 Theoretical description and mathematical problem. 

Consider the basic system of Euler equations for internal water waves in two dimensions (xz) with a stable 
stratified ambient state and the buoyancy frequency N(z). 

u x + w z = 

Po U t = -p (V, V) U - p x 

Po w t = PoJm, Y)w -p z - p g 

T' t +wT z =-(v,V)T' 

where u, v are velocity components, p a is the density, p is the pressure, p g is the body force due to 
stratifications, T z is the vertical background temperature gradient and T is the temperature variables. 

In the commonly accepted approximations (incompressibility, Boussinesq approximation, etc.) the 
solution of the system (|l[) is constructed as the following representation for the current function^ (x, z, t), 
(See, for details Appendix A) 

L 

i>(z,x,t) = Y,Z n (z)9 n (x,t), (2) 

n=l 

where Z n (z) are solutions of the correspondent Sturm-Liouville problem (Equation ( fL9|) in Appendix A) 

N 2 

Z zz + —Z = 0, Z(0) = Z(h) = 0, (3) 

and describe a vertical shape of the wave modes. The linear propagation velocities, c n , play the role of 
eigenvalues. The series (|2|) may be considered in the context of the Galerkin - Bubnov method, (for a 



(1) 



basis of this method see, for example, Gottlieb and Orszag (1989)). The coefficient functions 6 n (x,t) are 
solutions of the coupled KdV system (equation ( |2"7| ) in appendix A ) 



a? + cj n x + aJ2 a n m , k e m e k x + (3 2 d n e n xxx = o, (4) 

771, k 

a, f3 2 are scale parameters. Each solution 9 n (x, t) stands for n-th mode wave amplitude as a function of 
the space x G (—00, 00) and time t coordinates. Here the constants g m kn, d nm are nonlinear and dispersion 
constants. The coefficients g m kn, d nm have been derived by Leble (1991) to have the formulae (equation 
(EH) in appendix A ) 



h 

aN 2 c 2 

9m,k o 



1 2 \ _ 1 



— + - Z K Z™ z m z 

C m C k / c m c k 



z n dz, d n = ^- 2 (5) 



We specify the problem by zero boundary conditions for ip with respect to the variable z as follows from 
(§). Physically it corresponds to the zero value of the vertical velocity components at the upper (z = h) 
and lower boundaries (z = 0) levels. We select a localized along x axis initial condition described by a 
smooth enough function that model the paddle motion as in the experiments of McEwan (1983a). The 
function is also chosen antisymmetric along x in relation to the paddle axis centered in the middle of the 
tank. The time intervals of simulations are taken such that the initial disturbance decays essentially but 
does not reach the boundaries. 

In the simplest model to be presented here we consider the case of constant buoyancy frequency. In the 
concrete numerical calculations the choice of the constants is the following N = 1.23 s~ l for a stratified 
salted water in a rectangular tank 50 cm long filled to a depth 25 cm (Ref. McEwan (1983a)). 

4 Numerical method 

The numerical scheme is based on the solution of the nonlinear finite (difference) system 

((of 4 - (ni) i\ + c n - (ntx) /2h +Y.9 mk n (em ((o k y i+1 - mu) (6) 

k,m 

+ (d n - c n h 2 /Q) (>% 2 - 2 (ey i+1 + 2 {e-)U - (e n )j_ 2 ) / 2h 3 = o, 

where n, m, k are the modes numbers; i and j are discrete space and time variables respectively. The 
time step is denoted by r while the spatial step size by h. The equation (H) is accompanied with a discrete 
equation for the intermediate layer as: 



r|(0 n )'' +1 



')}) It + c n (V)ffi 1 - (0 n tf ) /2h + J29mkn (Of ' (V) - iP k t + }) l*h (7) 

k,m 



To support the results of further simulations to be presented here we would recover some mathematical 
ingredients of our work. Let {u n ){ = u n (xi,P) denotes the grid function obtained from the exact solution 
u n (x,t) of the Cauchy problem for the cKdV system calculated at the grid points Xi, P. Then the 



discrete solution (0 n )j of the system (]6|), (|7] ) defines the discrepancy vector (residual function) with 
the components (v n )i = (9 n )\ — («")-. Introduce the L>2 norm for V J ': 



The convergence of the scheme @, (0 ) at (r, /i — ► 0) follows from the 

Theorem 1 Let the solution u n (x,t) of the Cauchy problem for the cKdV system with t G [0,oo);x G 
(—00, 00) exists with the third continuous derivative with respect to time t, and the fifth derivative with 
respect to x. Then the following inequality takes place. 



&M\a n 

1-^/1 i^yW mo(t +h 2 ) 



V j+1 \\ < 1 = MO (r 2 + h 2 ) , (9) 
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where M is a constant independent ofr, h and #™ J = max ( <?^£ ); L is quantity of modes in (11). 



Therefore, the convergence follows from the inequality (||). The proof of the theorem is placed in the 
Appendix. We give the proof for more simple case then the only one of the numeric scheme equations 
(§) is used in simulations. The proof for the scheme (J§)-(0) with a half-step is analogous but significantly 
more cumbersome, and we decided to give not it here, taking in mind publishing in the forthcoming special 
mathematical paper . This scheme and proof may be of interest even in the scalar (single KdV) case in 
spite of existence of numerous versions of such results. We believe the ours is effective as it follows from 
the tests and experiments presented here. 

Remark 2 We have to note that evidently the existence of the fifth derivative by x and the third derivative 
by t are unnecessarily rigid conditions. These conditions are used for simplicity of consideration only, and 
may be easily weaken. One can replace them by the condition of the third order derivative by x and the 
first order derivative by x existence, then only the convergence is more slow. Moreover, one can accept 
a concept of weak solutions to be applied in applications. Usually such a concept allow to weaken the 
requirement of derivative existence on one order of derivatives; however the author have not verified it for 
the problem under consideration because necessity of such a delicate analysis is not evident. 

5 Computation and results 

When the coefficients in the Sturm-Liouville equation @ are constant it has very simple general solution 

Z n = B n Sm (2f ) , n=l,2,3,...,L 

which tends to zero at boundaries. The eigenvalues 

Nh , . 

c n = ,n = 1,2,3,. ..L. 10 

mi 

h 

have the sense of linear IGW velocities. Normalization is determined by J (Z n ) 2 N 2 dz = 1 and gives 





B n = (wh)' , Hence 



z"U)= I — \ 5 sin(^) in: 



h 



As it was discussed in the Sec. 2 one can select the initial perturbation for the current function (1) which 
have the general form at the point t = 

x, 0) = £ Z n (z)0"(:r, 0) = z) = ^i(x)^(z). 

n=l 

Taking smooth, exponentially localized initial condition for tpi(x) while ^(z) that model the impact 
of some rotating paddle motion as 

<fi(x) = Co g h (xj and ^(z) = (jfij^J Sech(b(z - z ))Tanh(b(z - z )) (12) 

The choice of the form and the constants a, b, 1 qualitatively reflects the paddle movement (we restrict 
the movement by some isolated pulse) and the numerical value for the amplitude is estimated also from 
the description of the experiment of McEwan (1983a). Then 

L 

(#,vo = 5^,innx,o) = (z j ,m^)) m*)- (13) 

n=l 
h 

[Z\ Z n ) = Jn 2 SP Z n dz = {1, (j = n)}, {0, (j ? n)} (14) 
o 

h 

(Zi,<p 2 (z)) = J N 2 Z^ 2 (z)dz. (15) 
o 

For h = 0.25m and N = 1.23s -1 , taking the first five nonzero modes with coefficients from ( |T5| ) and the 
numbers (2,4,6,8,10) while the odd numbers give zero, together with (|I4|) and putting in ([T3| ) to obtain 
9(x, 0). The estimation of the modes energy satisfactory keep that of the initial perturbation. Then using 
(1) we get the initial perturbation i[)(z,x, 0). We give cross section plot at the middle of the tank Fig. (2) 
to show the antisymmetry and the percentage of energy content in the initial perturbation model 




Fig.2 y?2(z) initially and after Fourier representation 



For the wave profile after time t we solve numerically system (|j) with coefficients from (||]), fllPf ) to obtain 
8 n (x,t) then using (1) for ip(z,x,t) evaluation. For the chosen constants the calculated coefficients of 
nonlinear, dispersion terms and modes propagation velocities are picked together in the table 
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d 22 =.00004, d 44 =.000005, d m =.000002, d 88 =. 000001, di io=-0000003, 
c 2 =.05, c 4 =.025 , c 6 =.016, c 8 =.012, ci =.0098 

For the case of linearized equations the wave profile after time t=0.02 s may be obtained directly from 
the same system if one put nonlinear constants to zero. The evaluation result is shown in the following 
contour plot Fig. (3) for the upper half of the tank. 
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Fig. 3 The wave profile 3-dimensional plot in linear case with dispersion for the upper half of the 

tank (12.5 cm in z and 50 cm in x directions) 



For the nonlinear case we show one- dimensional (x) profiles Fig.(4.a,b,c,d,e) for the second, forth, sixth, 
eighth and tenth modes respectively 
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Fig. (4. a) Second mode, initially (dotted) and at t=0.02 sec (continuous). 




Fig.(4.b) Fourth mode, initially (dotted) and at t=0.02 sec (continuous). 




Fig.(4.c) Sixth mode, initially (dotted) and at t=0.02 sec (continuous). 
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Fig.(4.d) Eighth mode, initially (dotted) and at t=0.02 sec (continuous). 




Fig.(4.e) Tenth mode, initially (dotted) and at t=0.02 sec(continuous). 
Fig. 4 One-dimensional (x) plots of the wave profile in nonlinear case. 



The contour plot for the wave profile at t — 0.02 sec is presented below (Fig. 5) for the same half of the 
tank. 
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Fig. 4 Three dimensional (x,y,z) plot of the wave profile for the upper half of the tank. 
The horizontal numbering indicates the number of points used in plot while the vertical 

show the actual z profile of the waves 

All the plots for the modes and for the sum (current function) show the behavior that is typical for the 
multisolitonic perturbation. It looks likely a decay of the initial condition to solitons in the single KdV 
equation theory. The difference between linear dispersion and nonlinear looks crucial. Even within the 
short time interval we present here the horizontal scale of the disturbance radically changes. The final 
scale is defined by constants of the cKdV system and the amplitude of initial condition. The process of 
the wave propagation is accompanied by interaction that imply the energy transfer between modes. This 
phenomenon may be considered as a possible reason for the vertical fine structure generation (Ref. Leble 
(1991) ) The combination of the fine structures may explain the experiments of McEwan (1983 a,b). 

6 Conclusion 

We conclude with the following remarks about restrictive assumptions we made and perspectives. First in 
the framework of clarity of the plotting and interpretation we considered only one direction of propagating 
waves. We would repeat that one can elongate the time of simulations to have more interactions but such 
picture looks cumbersome. 



The relatively small number of modes is taken into account because it conserves most the energy. The 
number may be increased but at the times we consider the effect of small- vertical-scales modes generation 
may be neglected. Returning to the general significance in the theory of the Euler equations validity 
mentioned in the introduction, one can easily put the important question. The question arises when we 
analyze the expressions for the constants in the of equations cKdV system. In a limit of the homogeneous 
medium, or N — > 0, the nonlinear constant increase while the dispersion and velocities go to zero. It 
means the growth of solitons number and more strong interaction between modes. The discussion of this 
phenomenon and ideas of adequate regularization will be considered in the forthcoming publication. 
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7 Appendix A 

Theoretical description and mathematical problem. 

Consider Euler equations that describe IGW in two dimensions (x,z) stratified water 

Ux + w z = 

PoU t = ~Po{v,V)U-px i , * 

po w t = p _(v, V)w-p z - p g 
T' t +wT z =-{v,V)T' 

where u, v are velocity components, p a is the density, p is the pressure, p g is the body force due to 
stratifications, T z is the vertical background temperature gradient and T is the temperature variables. 
Combining equations in (|Tj]) and using the state equation for liquid p = —p Q aT , a = — -=-is the 

coefficients of thermal expansion T , = — we obtain 

L * ag 



Aw tt + N 2 w xx - N 2 



(v,V) / wdt 



+ 



v, V) / w z dx 



+ [(v,V)w) txx = (17) 



txx 

txz 



Going to the dimensionless variables by the following formulas 

x t = \ lX >,t=2jL t > : u = w = i±JL w >, a V = T,N= |iV' 

So (dw/dxi) = (1/Ai) (dw/daQ , (dw/dt) = (W/3/2n) (dw/dt') 
Substitute in fli7| ) and omit primes for simplicity 



u 



Jq w x dt + w f* w z dt 







Introducing the current function , wt = —crijj x , u? = aip z , Integrate w.r.t. x 

ipzztt + N 2 tp xx = -(3 2 ip xxtt - cr [ip z ip xz - ip x ipzz) zt + <?N 2 



ipz I ^xxdt-ip x / ifj xz dt 
o Jo 



Substitute in by ^ = E z m m , multiply by z n , inegrate w.r.t. z and using The separation of 

m 

variables that gives 

N 2 

*« = -2T* n (19) 

C n 

So we have 



ell 
N 2 



nn _ 2 /in 



m.k 



Ki,k Q m ®t + e m,k 9™ I 8 K x dt 



h 



where a" (fc = iV 2 / + ^) z^zVz, 6" , fc = ^ / z™z fe z"tiz, C >fc = iV 2 / * 
-ft v fe m/ fe -ft -ft 

The equations of the separated propagated modes are obtained as follow 

Let 9? = u n , c n B n x =v n 

So (l20f) becomes 

c 3 /3 2 

m™ — c n w™ = -j^rv xxx + Nonlinear terras 



z m z^z n dz 



Using the projection operator P, P + 

p + ( u ' 

\ v 



v? - c n u x = 
1 



(20) 



(21) 



1 + £^,92 



1 _ s^ld 2 

1 2N 2 U x 



27V 2 u x 
1 



V = (ip n ' 



P_ 



if' 



P- = 

—kip' 

- ^-) 



_ 1 _ £^,92 

1 27V 2 u x 



.1 + ££1,92 
^ 27V 2 "a 



[22) 



Operating P + , P on (|2~ID and using ( p2|) we obtain the equations for the separated modes <p n+ , <p n 



- c n ^ x + - %^ x t x - ^e <, k j> + + dt 



m,k 



Cfe 



+ 



(ip+-ip ) m \.{-4>+~ip ) k 



'm,k 



I- 



dt 



(23) 



(f t -f C n (/9 X . -|- 2N?<P XXX 



m.k 



+ 



>" + -^-) m r (^ + -^ N ' 



/ 



ci- 



(24) 



Let 



So (E3 



becomes 



-\k 



+bl >k (*+ + vl/-) - + ^-^ + ^ (vp+ _ ($+ - vL 



(25) 



c 3 f3 2 
' 2N 2 



-6JL fc (*+ + vI/~) - 



Cfe 

m,fc 



(26) 



In next work we will consider both directed modes. In this work we evaluate only for short time to 
describe the phenomena, so we consider only one direction that have the form 
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8 Appendix B 

Proof of the convergence. 

Now we realize that the finite-difference method used gives us a correct solution of the problem if the time 
and space steps is enough small and we show how we can easily estimate an error of the numerical method. 
At first step we prove that a solution (6 1 ™)- 7 of finite-difference equations (H), |7|) converges to some solution 
of equation (??) as r, h — > 0. 

We give here a proof for more simple case for brevity when the scheme is based on equation @. The 
proof for equations (|j), [7]) is analogous but more cumbersome. Therefore, now we consider the difference 
equations 
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Here h is a space grid step, % is a discrete space variable, j is a discrete time, r is a time step, n is a mode 
number. We suppose that an exact solution is a smooth one and therefore 
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In this way, equations ( p9|) approximates equations (??) with an error of the order of 0(t + h 2 ). 

We know that a single KdV equation u t + uu x + u xxx = has the conservation law u 2 (x, t)dx = 
const. So we can consider a solution of a single KdV equation as an element of L 2 -space, parameterized 
by time t. Equations (??) are of the same structure as a simple KdV equation, so it is naturally to develop 
the proof in the functional space with L 2 -norm. 



Let (9 n ) 3 denote a column 
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and let Q 3 denote a column consisting of columns (9 n ) 3 introduced above: 

/ (9 l ) j \ 
(6 2 ) j 

& 3 = (9 3 ) 3 

{ W ) 

Let u n (x,t) be an exact solution of (??). Then (u™) 3 = u n (xi — ih,t — jr) is an exact solution in the 
point Xi of the grid at the moment jr. Let 
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The difference V 3 = Q 3 — U 3 is an error of the numerical method. The L2-norm of V 3 is 
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The difference solution (# n )^ converges to the exact solution (u 71 )^ if || V^\\ — > as r, h — > 0. 
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We substitute (??) into (|29| ) and obtain an equation for the error (v 
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Pick out a linear part of expression ([32]) and introduce for convenience an operator T J by the expression 
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Because (« n )^ is an exact solution of differential equations, the right part of ( j3~2"l) is a small of the order 
O (r + /i 2 ) . Therefore, we can rewrite (|32|) in the following way 
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We need an estimate of a norm of TK It is valid 
Lemma 3 If u n (x,t) is a function with a bounded derivative \-^u n (x, t)\, then for any b > and for all 

r<bh 



(34) 



there exists a constant A, independent of r, h, j, such that ||T J '|| S ^ exp(Ar). 



Here HT^'U is a spectral norm of the operator T- 7 , i.e., the norm of T J that is induced by L 2 -norm 
|| V^'H- We will omit the subscribe "s" in the following and will simply write \\T^\\. The proof of Lemma |3] 
is a standard one, but cumbersome. Therefore we do not give the proof of Lemma |3] here and recommend 
reader books [0, |f2T . 



Let F J be the column constructed from (/") J ' in the same manner as J constructed from (6™Y . We 
need a norm of F- 7 ; simple estimation gives 
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So, using the operator T J , we can rewrite fl32|) as 



The following estimate is valid 
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If we consequently substitute (^) into itself, we get 
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We have used a known formula for the sum of a geometric series. Here t max is the time of simulation: 

^ t ^ ^max • 

Substitution of instead of Hi^H into (pT7|) yields the inequality 
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Because of the initial conditions we know exactly, we have (t>™)° = 0. So, taking into consideration 
\V°\\ = we obtain the solution of the inequality 
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i.e. H^^ 1 !! — > as r, fa — > 0. The convergence is proved. 

The error of the method considered is of order O (r + fa 2 ), if the solution under consideration is enough 
smooth one. The method converges for any r, fa — * 0, but the rate of converging is dependent of relation 
between r and fa: the more p- the better. If the solution is not smooth, then we must consider the solution 
in the sense of a distribution theory; this case is not study yet. 

The proof of the convergence for the scheme (|6]), (|7| ) is analogous, but more cumbersome and tedious. 
The error of the method (|), (|7| ) is O (r 2 + fa 2 ), and instead of condition ([34] ) the condition r < bh A is 
used. The main idea of the proof is the same one: at small r, fa the error due to difference approximation 
of the nonlinear terms influence on the result more less then the errors due to difference approximation of 
dispersion terms; and all difficult estimates only are necessary to show this fact. Therefore, choice of the 
optimal step tau may be accomplished with the help of consideration of simple difference equation with 
dispersion term only. We had used this fact when made numerical experiments. 



